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The inspiral of two compact objects in gravitational wave astronomy is described by a post- 
Newtonian expansion in powers of (v/c). In most cases, it is believed that the post-Newtonian 
expansion is asymptotically divergent. A standard technique for accelerating the convergence of a 
power series is to re-sum the series by means of a rational polynomial called a Pade approximation. 
If we liken this approximation to a matrix, the best convergence is achieved by staying close to a 
diagonal Pade approximation. This broadly presents two subsets of the approximation : a super- 
diagonal approximation P^f and a sub-diagonal approximation where M = N + e, and e 
takes the values of or 1. Left as rational polynomials, the coefficients in both the numerator 
■ and denominator need to be re-calculated as the order of the initial power series approximation 
is increased. However, the sub-diagonal Pade approximant is computationally advantageous as it 
, can be expressed in terms of a Gauss-like continued fraction. Once in this form, each coefficient in 
the continued fraction is uniquely determined at each order. This means that as we increase the 
order of approximation of the original power series, we now have only one new additional coefficient 
to calculate in the continued fraction. While it is possible to provide explicit expressions for the 
continued fraction coefficients, they rapidly become unwieldy at high orders of approximation. It 
' is also possible to numerically calculate the coefficients by means of ratios of Hankel determinants. 

However, these determinants can be ill-conditioned and lead to numerical instabilities. In this 
article, we present a method for calculating the continued fraction coefficients at arbitrary orders of 
^q' approximation. 

I. INTRODUCTION. 

> 

, Gravitational wave (GW) astronomy is one of the newest and most exiting fields in astronomy and astrophysics. 
At present, a number of the ground based interferometers, such as Advanced LIGO and Advanced Virgo, are being 
upgraded with the goal of being back online in 2015 [l|, Q. There are also plans for a cryogenic detector in Japan 
^vq . called LCGT [3] and a large subterranean detector in Europe called The Einstein Telescope |4j. In order to detect 
t— ( ' more massive sources of GWs, there are also plans for an ESA led spaced based GW mission called eLISA [5|. One 
of the most powerful sources of interest is the inspiral and merger of compact binary systems such as binaries of 
supermassive black holes from galaxy mergers, or the inspiral of stellar mass black holes and neutron stars. In all of 
these cases, the rate of energy loss due to the emission of GWs is given by the energy balance equation 



X 



dE(v) 

m—jj— = -F(v), (1) 

where m is the total mass of the binary system, E(v) is the binding energy of the system, v is the linear velocity 
of the system components and F(v) is the GW flux. Most of the search methods in the field are based on optimal 
Wiener or matched filtering, where phase matching between the model template and a possible source in the data is 
the highest priority. As the phase of the GW is a function of both E(v) and F(v), in order to successfully search for 
these systems in a data set, we need to be able to successfully model the GW flux function. However, as we do not 
have an exact solution for the GW flux from a binary system, it is represented by a post-Newtonian expansion (i.e. 
an expansion in powers of (v/c), where c is the speed of light). This series has the form 



F(v) = Fj 



N 



(2) 



.4=0 i=6 

where n has different values depending on the type of physical system in question. The quantity Fn = 32ij 2 v 10 /5 is 
the Newtonian flux and rj — miTO2/m 2 is the symmetric mass ratio. As the post-Newtonian approximation is based 
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on the limiting condition v <C c, it is thought that the above representation of the flux is asymptotically divergent as 
we approach the domain where v ~ c. An improved model, based on Pade approximation, was proposed where the 
flux is written in the form I™ 



F p (v) =F N (1- v/vi r y 



\n(v/v lso )y^JiV l 



l_i=0 



(3) 



where the coefficients are defined by ij = li(Ai, B{) and a, = di(Ai, B il vi so ). This flux representation was demonstrated 
to have better convergence properties than the PN flux for both spinning and non-spinning test mass systems @, @] , 
and also for non-spinning comparable mass binaries Q. In Equation the quantities (vi so ,vi r ) correspond to the 
orbital velocities at the last stable and unstable circular orbits. 

In approximation theory one can equate a Pade approximation to a power series by 



JV 

Vrj v - p n = 1=0 



(4) 



i=0 



such that the n + 1 coefficients on the left hand side are matched by the N + M + 1 coefficients on the right hand side. 
Theory suggests that the best convergence is obtained when the Pade approximation is kept to near diagonal values, 
i.e. N » M. This allows us to define two types of approximant : sub-diagonal Pade approximants (M > N) and 
super-diagnoal Pade approximants (N > M). The advantage of defining the two types of approximants is that we 
sometimes obtain poles in the Pade approximation. However, a pole in say the sub-diagnoal approximant is usually 
not reflected in the super-diagonal approximant. In the particular case where M = N or M = N + 1, we can write the 
Pade approximant in a continued fraction form. The advantage of this is in the computation of coefficients. In the 
case of a super-diagonal approximant, each time we increase the order of the power series, we have to re-compute all 
N + M + 1 coefficients in the Pade approximation. However, in the sub-diagonal case, where we have re-written the 
approximation as a continued fraction, at each new order in the power series, we only have to calculate the newest n th 
order coefficient, as all others stay constant. We can now define the operator C n [...} in Equation Q as the continued 
fraction representation of a sub-diagonal Pade approximant, i.e. 



i=0 



CO 



(5) 



c 2 x 



1 + c n x 



The continued fraction used for the flux model corresponds to a sub-diagonal Pade approximant with M = N + e 
and e having values of or 1, i.e. 

JV 



On 



i=0 



p 



N _ i=0 



(6) 



While the coefficients in the continued fraction can be solved for in terms of explicit expressions using a computer 
algebra system such as MAPLE or Mathematica, it is difficult to present the expressions beyond a certain order due 
to their length and complexity. Our goal here is to present a fast numerical algorithm to calculate the coefficients of 
a continued fraction at arbitrary order n. 



II. THEORETICAL FOUNDATIONS 



In cases where a power series has slow convergence, a rational function approximant such as Pade approximation, 
and by extension a continued fraction representation, have been shown to increase the speed of convergence @. In 
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general, a power series converges within a tight radius of convergence, but only close to the point of expansion and 
in between the first poles. A power series can also be numerically unstable due to the fact that, as we increase the 
number of terms in the expansion, the magnitudes of the exponents of the monomials continue to grow. A continued 
fraction, on the other hand, converges everywhere in the complex plane, except at the poles. However, as the continued 
fraction a ppr oximation to a power series contains a monomial in each denominator, they are quite adept at capturing 
the poles [loj]. 

While continued fractions have a number of attractive properties, there are certain numerical issues in calculating 
the continued fraction coefficients, given a generatin g p ower series. We will briefly outline this in the following 
example. Firstly, by using the following nomenclature [llf 



i=l Oi 



ai 



<i 2 



03 



(7) 



we can define a regular continued fraction (from here-on, C-fraction) [jj 



oo „. r a(i) 



(8) 



where a(i) is a constant. If we now take a(i) = 1, and equate the C-fraction with a power series, we can write 



(9) 



In this case, the coefficients of the C-fraction, ai, are given by [lj 



(i) 



«2m — 



tt(1) rr(2) 
n m-\ n ™ 

11 ™ n m-\ 



1 0-2m+l 



rrU) rr(2) 
JJ m+l J7 m-l 

o-(l)o-(2) ' 



(10) 

(ii) 



where we define the Hankel determinant 



HW=HW(a) = det 



^ Cn ' * ' Cm-\-n—2 Cm-\-n — l \ 

Cn+1 C-n-\-2 ' ' ' ^m-\-n—l Cm-\-n 

\Cm-\-n— 1 C7n-\-n ' ' ' C2rn-\-n — 3 C2m+n — 2 J 



(12) 



The problem with this method of finding the C-fraction coefficients is that the Hankel determinants are frequently 
ill-conditioned, leading to a number of numerical issues. We can see from the above expression that in the case of one 
of the Hankel determinants in the denominator being very small, the C-fraction coefficient blows up. 



III. CALCULATING THE RATIONAL FUNCTION COEFFICIENTS. 



Our goal in this work is to present a method that avoids using ratios of determinants and will provide a more 
numerically stable method of calculating the C-fraction coefficients. The first step in the analysis is to have expressions 
for the coefficients of the numerator and denominator in the rational polynomial defined by Equation (JB]) at arbitrary 
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order. By expanding at lower orders of n, we find that the first five rational function expansions have the form 
Co 



PI' 



P? 



1 + C\X ' 

cO (1 + c 2 x) 
1 + (ci + c 2 )x'' 

+ c 2 + r 

cO (1 + (c 2 + c 3 + c 4 ) x + C2C4X 2 ) 



A 1 = 

J- T 

pi = = ;°( 1 + ^+ c 3)^) (13) 

f + (ci +c 2 + c 3 )a; + ciC3a; 2 



f + (ci + c 2 + c 3 + c 4 ) x + (ci (03 + c 4 ) + c 2 c 4 ) a; 2 ' 

cO (l + (c 2 + c 3 + c 4 + c s ) x + (c 2 (c 4 + c 5 ) + C3C5) x 2 ) 



f + (ci + c 2 + c 3 + c 4 + c 5 ) x + (ci (c 3 + c 4 + c 5 ) + c 2 (c 4 + c 5 ) + C3C5) a; 2 + cic 3 c 5 a; 3 ' 
By investigation, the di coefficients in the numerator follow the form 

da = c , 

AT+M n 

di = c E Cj = co^Cj, 

n — 2 n 

rf 2 = Q^E E C,'Cfc, (14) 
j=2 fc=j+2 
n — 4 n — 2 n 

rf 3 = c °E X! X! c i c * c i> 

j=2 k=j+2 l=k+2 

where each di contains i internal summations. By generalising to arbitrary order, we can express the coefficients in 
the numerator of the rational function in the form 

h ri+2 n+4 n— 2 n 

4=C0^ E E •■• E X! c J c fe c i--- c rC s , (15) 
j=2 k=j+2 l=k+2 r=q+2 s=r+2 

where n = n — 2{i — 1). Each outer summation begins with a lower limit of j = 2 and an upper limit of n — 2(i — 1), 
while each inner summation has a lower limit of s = r + 2 and an upper limit of n, where the summation index r has 
an initial value of r = 2i. 

Similarly, by investigation, the coefficients in the denominator of Equation © follow a similar pattern where 

fa = 1, 

N+M n 

fi = E c j =E C ." 

3=1 3=1 

n — 2 n 

h = E E c ^> ( 16 ) 

j=l fe=j+2 
ro— 4 7j— 2 n 

/ 3 = E E E c J cfeQ - 

3=1 fc=j+2 Z=fc+2 

(17) 

where for each /j, there are i internal summations. Once again, generalising to arbitrary order, we can write the 
coefficients of the denominator in the form 

n n+2 n+4 n — 2 n 

fi = E E E - E E c J c kci-c r c s . (is) 

j=l fc=j+2 Z=fc+2 r=q+2 s=r+2 

where again n = n — 2(i— 1). This time the lower limit of each external summation begins at unity, with an upper limit 
again of n — 2(i — 1), while the internal summation again has lower and upper limits of s = r + 2 and n respectively, 
where in this case, the summation index r has an initial value of r = 2i — 1. 
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IV. CALCULATING THE C-FRACTION COEFFICIENTS. 



Now that we have the coefficients for the rational function, we can equate the original power series with the rational 
function, i.e. 



N I M 



i=0 i=0 I i=0 

By multiplying both sides by the denominator on the right hand side, we define a new equality 

n N 



J2hx l =J2d t x\ (20) 



i=0 i=0 

where the series on the left hand side is a truncated series of order n, i.e. 

n n M 



^&^ = ]T a ^]r/ 4 x 4 +0(x" +i ), (21) 

i=0 i=0 i=0 

and where we define the coefficients 6, by 

v 

h = a t + ^2ai-jfj. (22) 
i=i 

Here y = min(i, \n/2~\), where the ceiling \x~] of any real number x is the unique integer m satisfying \x~] = m <s=>- 
x < m < x + 1. Thus for each order of approximation n we obtain a system of n equations for each bi of the form 
hi = di. For i — we have the trivial solution b = do = c . However, at higher orders, there are a number of 
subtleties that allow us to simplify the above system of equations. 

If we look at the case of b\, when n = 1, we end up with the expression a\ + ao/i = ai + aoci = d\ = 0, with the 
simple solution c\ = — ai/ao. Now, if n = 2 we obtain a similar expression, but in this case d\ is non-zero, i.e. 

2 2 

ai + ao/i = ai + a c,- = c c,-. (23) 

Using the fact that ao = Co, this expression again reduces to ai + a ci = 0. In general, as we go to higher orders of 
n, we end up with an expression of the form 

n n 

bi = a x + a Q ^ c 3 = c o X! C J 

i=i i=2 

=> ai + a ci =0 V n > 2. (24) 

As a further example, we look at the coefficient bi. In this case ifiV<i = 2, rfj = 0. Thus at n = 2 and n = 3, we 
end up with the equations 

a 2 + ai/i = 0, (25) 

a 2 +a 1 f 1 +a Q f 2 = 0. (26) 

The first equality for n = 2 provides the solution a 2 + a\{ci + c 2 ) = 0, while for n = 3, we obtain 

3 13 

a 2 + ai/i + ao/2 = «2 + Qi ^ Cj + a CjCfe = 0, 

j = l j = \ k=3 

3 

=> a 2 + ai ^ Cj + a cic 3 = a 2 + ai(ci + c 2 ) = 0, (27) 
i=i 
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as aoCi — — a\. At n = 4, we obtain our first non-zero value of d,2 giving the equation 02+01/1+00/2 = e?2- Expressing 
both sides in terms of Oj's and Cj's we get 



4 2 4 2 4 



02 + ai 2jcj + a ^ ^ CjCfc = cp 2J c J c fc- ( 28 ) 

3=1 3 = 1 fe=3+2 3=2 fe=j'+2 



Subtracting the quantity on the right hand side gives us the expression 

4 1 4 

g-2 + qi y^ Cj + ap y^ c 3 c fe = o 

j=l 3=1 fc=J+2 

4 4 

=> a 2 + ai 2J Cj + a ci ^ ct = a 2 + a\{ci + c 2 ) = 0, (29) 

3 = 1 fe=3 

where again we have used the fact that aoCi = —a\. Thus, we can see that regardless of the value of n, we end up 
with the same equation for 61 . It is also clear from Equation (f2"2"j) that initially each 6, is a function of the coefficients 
Cj, V < j < n through the functions /j. However, as we have demonstrated in the examples shown above (and 
others not displayed here), cancellations between the Cj's manage to kill the higher order coefficients at i < j < n, 
meaning that finally, the bi coefficients are a function of the coefficients Cj, V < j < i. Thus, we can now write the 
system of equations at each order in the form 

bt = ai +J2 ai-jfj = 0, (30) 
3=1 

where 

i i+2 i+4 i-2 i 

ij = ^2 12 12 ■■■ 12 12 c kcic m ...c r c s , (31) 

fe=l l=k+2 m=l+2 r=q+2s=r+2 

and i = i — 2(J — 1). As this system of equations are linear in Cj, we can express Equation (|30[) in the form 

m + fad = 0. (32) 

By factoring out the ct coefficient in Equation (1221) and by equating this expression to zero, we can define a system 
of iterative linear equations 

\i/2] li/2\ / L(i"l)/2J \ 

a,i + ^2 a i-jfj = a i + 12 ai -ofj + c i a (i-i) + 12 a (i-i)-jfj = V i > 1. (33) 

3 = 1 3 = 1 \ 3 = 1 / 

Here the floor lx\ of any real number x is the unique integer satisfying [x\ = m x — 1 < m < x. In the above 

expression we define the reduced coefficients fj as 

i i+2 1+4 (i-l)-2 i-l 

fi = 12 12 12 - 12 12 c kcic m ...c r c s (34) 

k=l l=k+2 m=l+2 r=q+2 s=r+2 

where i = (i — 1) — 2(j — 1). While these coefficients have the same form as the /i coefficients given in Equation (|18j). 
we should note that the upper limit is now a function of (i — 1) instead of n. As an example, the /, coefficients to 
order n = 3 are given by 

i-l (»-l)-2 (i-l) (i-l)-4(i-l)-2 (i-l) 

fo = 1 , fi = 51 Cfe ' f 2 = 12 X! CkCl ' ^ 3 = X! X] X! c ^ c ' c ™ ' ete - ( 35 ) 

k=l fe=l Z=fc+2 fc=l i=fc+2 m=i+2 



We also define the double- reduced coefficients fj, with again the same form as the coefficients in Equation (|18|) . except 
this time the maximum upper limit is given by (i — 2) instead of n, i.e. 

i i+2 i+4 (i-2)-2 i-2 

fo = 12 12 12-12 12 c kcici...c r c s (36) 

k=l l=k+2 m=l+2 r=q+2 s=r+2 
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where i = (i — 2) — 2(j — 1). We can see that the functions fj has the same form as the functions fj where we replace 
(i — 1) by (i — 2). Again, to n = 3 order, the coefficients are described by 

i-2 (*-2)-2 (i-2) (»-2)-4(*-2)-2 (i-2) 

/o = i , h = ^2 ck > = X! X! cfcQ > / 3 = X X! c t c ' c ™ ' etc - ( 37 ) 

fc=l fe=l l = k+2 fe=l (=fc+2 m=i+2 

Finally, solving for Cj, we can define the expression 

L«/2J 

Of + "■ 

Cl = = LC^Iyi] / (38) 

3=1 

We can see from the above expression, that the numerator and denominator of the a coefficient have a similar form, 
but with the index i in the numerator, and (i — 1) in the denominator. It is not difficult to then show that the 
denominator of the Ci coefficient is identical to the numerator in the Cj_x coefficient. Therefore, at the i th order we 
only need to calculate the numerator at that order, i.e. 

[»/2J 

OLi = a,; + ^2 a i~]fji (39) 
3=1 

and can thus define 

Ci = — (40) 

Q!t-1 

as the coefficients of a C-fraction representation of a sub-diagonal Pade approximation at arbitrary order. We have 
tested the accuracy of the method against value obtained from the analytic expressions and have found perfect 
agreement. 



V. CONCLUSION 



In this work we have looked at a way of calculating the coefficients of a C-fraction representation of a sub-diagonal 
Pade approximation at arbitrary order. Pade approximation has become a common method in improving the conver- 
gence of a post-Newtonian series for binary compact objects. A sub-set of the Pade approximant, called a sub-diagonal 
Pade, is easily expressed as a regular continued fraction, called a C-fraction. The advantage of the C-fraction repre- 
sentation is that, as we expand the generating power series to higher orders of approximation, we only have one new 
coefficient to calculate. While it is possible to analytically derive expressions for the coefficients of the C-fraction, 
these expressions become very lengthy even at low orders of approximation. Numerical methods exist that are based 
on calculating the ratios of Hankel determinants. However, if the determinants are ill-conditioned, it makes the ap- 
proach unstable. We present here an iterative numerical method that, given the coefficients of the generating power 
series, allows one to calculate the C-fraction coefficients at arbitrary order 
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